library(ncdf4)
library(ggplot2)
require(data.table)

#change to full local path to the replication data
wd1 <- "C:\\Users\\YFan\\OneDrive - NORCE\\NFood_replication\\"

#=================estimate effect of volcanic eruption on yield using SAOD from El Chichón and Pinatubo eruptions

#according to Proctor et al., Pinatubo eruption (global average of +0.15 SAOD) reduced direct sunlight by 21%, 
#increased diffuse sunlight by 20% and reduced total sunlight by 2.5% (Fig. 1f,  Extended Data Table 1)
#verify by translating log effect back to % using coefficients in their Extended Data Table 1
dAOD = +0.15; #after Pinatubo
(exp(-1.406*dAOD)-1)*100; (exp(1.171*dAOD)-1)*100 #-19% direct / +19.2% diffuse if using AOD coefficient from Pinatubo
dAOD = +0.09; #after El Chichón
(exp(-1.039*dAOD)-1)*100; (exp(2.063*dAOD)-1)*100 #-9% direct / +20.4% diffuse if using AOD coefficient from El Chichón 

dAOD = 0.084 # for the SAI - RCP45 stablizing experiment in Proctor
(exp(-1.406*dAOD)-1)*100; (exp(1.171*dAOD)-1)*100 #-11% direct / +10.3% diffuse if using AOD coefficient from Pinatubo 
(exp(-1.039*dAOD)-1)*100; (exp(2.063*dAOD)-1)*100 #-8% direct / +19% diffuse if using AOD coefficient from El Chichón 


#==============check the AOD dynamics of NorESM for historical period (1980-2014)
#TOM had all the variables I need and the top of the model is very close to the top of the stratosphere, this should be OK for my purposes.  
#Also, the "Description of the NCAR Community Atmosphere Model (CAM 3.0) ", 2004, states on p. 113 that "…solar fluxes are calculated… for the top of model (TOM) at layer 1 and the top of atmosphere (TOA) corresponding to layer 0.  
#The TOM fluxes are used to compute the model energetic balance, and the TOA fluxes are output for diagnostic comparison against satellite measurements."
#SOLIN=FSDTOA=FSNTOA+FSUTOA
#TOM radiative imbalance: FSNT(shortwave)+FLNT(longwave)
#Diffuse radiation at surface: SOLSD(vis)+SOLLD(nir) or FSDSVI+FSDSNI
#SAOD: TAUA550 or TAUA550

#========AOD variables from NorESM1-ME:
#TAUA550 "Aerosol absorptive optical depth at 550nm"
#TAUE550 "Aerosol optical depth at 550nm"
#or NorESM2: 
#AOD_VIS "Aerosol optical depth at 0.442-0.625um", 
#AODVVOLC "CMIP6 volcanic aerosol optical depth at 0.442-0.625um" (specifically for volcanos)
#these AOD variables need to be scaled by daylight fraction

solar.hist <- nc_open(paste0(wd1, "/data/climate_data/NorESM1-ME-HIST.2deg.cam2.h0.1979-2005-solar.nc"))
syr = 1979
fyr = 2005
nyr = fyr-syr+1
nlat = 96
nlon = 144

names(solar.hist$var)
names(solar.hist$dim) #dimensions "time" "cft"  "lat"  "lon"
dayfoc <- ncvar_get(solar.hist, "DAYFOC") #daylight fraction to rescale AOD 
aod.vis <- ncvar_get(solar.hist, "TAUE550") #AOD for visible band (550nm)
aod.abs <- ncvar_get(solar.hist, "TAUA550") #absorptive AOD
summary(aod.vis);summary(aod.abs) #aod.vis is the correct variable to use!

#fsdtoa <- ncvar_get(solar.hist, "FSDTOA") #Downwelling solar flux at top of atmosphere (FSDTOA=SOLIN)
#fsntoa <- ncvar_get(solar.hist, "FSNTOA") #Net solar flux at top of atmosphere (FSNTOA=FSDTOA-FSUTOA)
#solin <- ncvar_get(solar.hist, "SOLIN") #total solar insolation at TOA (=FSDTOA, but /= FSDS in the lnd model)
fsds <- ncvar_get(solar.hist, "FSDS")     #Downwelling solar flux at surface (even lower than FSNTOA)
solsd <- ncvar_get(solar.hist, "SOLSD") #Solar downward visible diffuse to surface (SOLSD=FSDSVI for vis band, SOLLD=FSDSNI for nir)
solld <- ncvar_get(solar.hist, "SOLLD")
summary(fsds);summary(solsd);summary(solld)
nc_close(solar.hist)

gridmx <- array(1:(nlat*nlon),varsize[1:2]) #matrix of grid index
solar.mon <- data.table(aod.vis=c(aod.vis/(dayfoc+1e-10)), fsds=c(fsds), solsd=c(solsd), solld=c(solld),
                        Sample=rep(gridmx, time=nyr), #grid indices
                        Year=rep(rep(syr:fyr,each=12),each=nlat*nlon), Month=rep(rep(1:12,time=nyr),each=nlat*nlon))
solar.m <- solar.mon[,.(aod.vis.m=mean(aod.vis),fsds.m=mean(fsds),solsd.m=mean(solsd), solld.m=mean(solld)),
                     by=.(Month,Sample)]
solar.dm <- solar.mon[solar.m, .(Year, aod.vis, fsds, solsd, solld,
                                 aod.vis.dm = aod.vis-i.aod.vis.m),
                      by=.EACHI, on=.(Month,Sample)] #de-mean by subtracting the mean of all observations from the respective grid and month-of-year)

solar.yr <- solar.dm[,.(aod.vis=mean(aod.vis),fsds=mean(fsds),solsd=mean(solsd), solld=mean(solld),
                        aod.vis.dm=mean(aod.vis.dm, na.rm=T)),  by=.(Year,Sample)]
solar.gl <- solar.yr[,.(aod.vis=mean(aod.vis),fsds=mean(fsds),solsd=mean(solsd), solld=mean(solld),
                        aod.vis.dm=mean(aod.vis.dm, na.rm=T)),
                     by=.(Year)] #this is unweighted by crop area!

plot(solar.gl$solsd~solar.gl$Year, type="l") #El Chichón (1982.3) and Pinatubo (1991.6) eruptions is evident from diffuse radiation change
plot(solar.gl$solld~solar.gl$Year, type="l") #but not evident in SOLIN (does not equal to FSDS in the lnd model), should use FSDSVD+FSDSVI from lnd model
plot(solar.gl$fsds~solar.gl$Year, type="l")  #FSDS (surface) /= SOLIN (TOA)
plot(solar.gl$aod.vis~solar.gl$Year, type="l") #AOD peaks after 1982 and 1991 very similar to Protocor, but values in other years are higher than 0
plot(solar.gl$aod.vis.dm~solar.gl$Year, type="l")

summary(solar.gl$aod.vis) #mean ~0.13 (both tropospheric and stratospheric AOD)
solar.gl[Year==1992, aod.vis] #0.217
solar.gl[Year==1983, aod.vis] #0.177 same as aod.vis in NorESM2
solar.gl[,.(Year, dAOD=aod.vis.dm)] #deviation the mean 0.085 for 1992, 0.045 for 1983 (similar to aod.volc of NorESM2)
#dAOD is a good estimate of volcanic only AOD, similar to aod.volc from NorESM2



#=========use NorESM2 (CMIP6) data========
syr = 1979
fyr = 2009
#time period of Proctor study (1979-2009)
nyr = fyr-syr+1
nlat = 192
nlon = 288
solar.hist2 <- nc_open(paste0(wd1, "/data/climate_data/NorESM2-MM-HIST.1deg.cam.h0.1975-2014-solar2.nc"))
aod.hist <- nc_open(paste0(wd1, "/data/climate_data/NorESM2-MM-HIST.1deg.cam.h1.1979-2009-AOD.nc")) #annual data, already rescaled by daylight fraction

names(solar.hist2$var)
names(solar.hist2$dim) #dimensions "time" "cft"  "lat"  "lon" 
v1 <- solar.hist2$var[[1]] #read all info of a variable
v1$name;v1$units;v1$longname #time is always the last dimension
solar.hist2$var[[2]]$longname
(varsize <- v1$varsize)

dayfoc <- ncvar_get(solar.hist2, "DAYFOC", start=c(1,1,(syr-1975)*12+1), count=c(nlon,nlat,nyr*12)) #Daylight fraction (used to rescale AOD)
aod.vis2 <- ncvar_get(solar.hist2, "AOD_VIS", start=c(1,1,(syr-1975)*12+1), count=c(nlon,nlat,nyr*12)) #AOD for visible band (0.442-0.625um)
aod.volc <- ncvar_get(solar.hist2, "AODVVOLC", start=c(1,1,(syr-1975)*12+1), count=c(nlon,nlat,nyr*12)) #CMIP6 volcanic aerosol optical depth at 0.442-0.625um"
# fsds <- ncvar_get(solar.hist2, "FSDS", start=c(1,1,(syr-1975)*12+1), count=c(nlon,nlat,nyr*12))     
# solin <- ncvar_get(solar.hist2, "SOLIN", start=c(1,1,(syr-1975)*12+1), count=c(nlon,nlat,nyr*12))     
summary(aod.vis2) #similar mean (~0.074) and 1st/3rd Qu. as NorESM1
summary(aod.volc) #mean ~0.0084 (max 024), much smaller than aod.vis2
nc_close(solar.hist2)

v1 <- aod.hist$var[[1]] #read all info of a variable
v1$name;v1$units;v1$longname;v1$varsize #time is always the last dimension
aod.vis1 <- ncvar_get(aod.hist, "AOD_VIS1", start=c(1,1,1), count=c(nlon,nlat,nyr)) #rescaled by daily data of DAYFOC
nc_close(aod.hist)

gridmx <- array(1:(288*192),varsize[1:2]) #matrix of grid index
solar.mon <- data.table(aod.vis1=c(aod.vis1[,,rep(1:nyr, each=12)]), #annual to monthly
                        aod.vis2=c(aod.vis2/(dayfoc+1e-10)), aod.volc=c(aod.volc/(dayfoc+1e-10)), #fsds=c(fsds), solin=c(solin), 
                        dayfoc=c(dayfoc), Sample=rep(gridmx, time=nyr), #grid indices
                        Year=rep(rep(syr:fyr,each=12),each=nlat*nlon), Month=rep(rep(1:12,time=nyr),each=nlat*nlon))
# solar.m <- solar.mon[,.(aod.vis.m=mean(aod.vis),aod.volc.m=mean(aod.volc),fsds.m=mean(fsds),solin.m=mean(solin)), 
#                      by=.(Month,Sample)]
# solar.dm <- solar.mon[solar.m, .(Year, aod.vis, aod.volc, fsds, solin, aod.vis.dm = aod.vis-i.aod.vis.m), 
#                       by=.EACHI, on=.(Month,Sample)] #de-mean by subtracting the mean of all observations from the respective grid and month-of-year)
solar.yr <- solar.mon[,.(aod.vis1=mean(aod.vis1), aod.vis2=mean(aod.vis2), aod.volc=mean(aod.volc), #no daily data of AODVVOLC for rescaling by daily DAYFOC
                         dayfoc=mean(dayfoc)),  by=.(Year,Sample)]

solar.gl2 <- solar.yr[,.(aod.vis1=mean(aod.vis1), aod.vis2=mean(aod.vis2), aod.volc=mean(aod.volc),
                         dayfoc=mean(dayfoc)), by=.(Year)] #this is unweighted and unmasked by crop area!

#plot(solar.gl2$solin~solar.gl2$Year, type="l") #but not evident in SOLIN (does not equal to FSDS in the lnd model), should use FSDSVD+FSDSVI from lnd model
#plot(solar.gl2$fsds~solar.gl2$Year, type="l")  #FSDS (surface) /= SOLIN (TOA) 
plot(solar.gl2$aod.vis1~solar.gl2$Year, type="l") #rescaled by daily data of AOD_VIS and DAYFOC 
plot(solar.gl2$aod.vis2~solar.gl2$Year, type="l") #AOD peaks after 1982 and 1991, but values in other years are high ~0.14
plot(solar.gl2$aod.volc~solar.gl2$Year, type="l") #0.054 in 1992

summary(solar.mon$dayfoc)  #monthly mean daylight fraction
summary(solar.gl2$dayfoc)  #mean daylight fraction 0.5
solar.gl2[Year==1992, aod.vis2] #0.224
solar.gl2[Year==1983, aod.vis2] #0.177
solar.gl2[Year==1992, aod.volc] #0.097 (<0.15)
solar.gl2[Year==1983, aod.volc] #0.052 same as noresm1-me
solar.gl2[,.(Year, aod.volc)]  

nc_close(solar.hist2)
rm(solar.mon, solar.yr, solar.hist2, fsds, solin)


#====all data masked and weighted by crop area below
c.hist <- nc_open(paste0(wd1, "/data/climate_data/NorESM2-MM-HIST.1deg.clm5.h0.1975-2014.nc"))
y.hist <- nc_open(paste0(wd1, "/data/land_cover/NorESM2-MM-HIST_crop_data.1975-2014.nc")) 
#no crop specific yield data available for the 8 crop types (all 16 cfts were simulated as generic C3/C4 crop in NorESM2-MM-HIST); 
#and not possible to separate climate/aerosol effects from land use effects in the historical simulation due to transient land cover and nitrogen fertilizaton rates!
area <- ncvar_get(y.hist, "ALL_AREA", start=c(1,1,syr-1975+1), count=c(nlon,nlat,nyr)) #annual data [lat,lon,years]
cft.area <- ncvar_get(y.hist, "CFT_AREA", start=c(1,1,1,syr-1975+1), count=c(nlon,nlat,16,nyr)) #16cfts

samples <- area>=0.001 #Mha -> 1000ha
samples[is.na(samples)] <- FALSE
samplesize <- apply(samples, 3, function(x)sum(array(1, list(nlon,nlat))[x], na.rm = T)) #different number of samples per year
samplesize #around 8100~8200 grid cells with generic crop
#masked by all crop area per grid
hist.y <- data.table(Year=rep(syr:fyr,each=288*192)[samples], 
                     area=c(area[samples]),  
                     Sample=rep(gridmx, time=nyr)[samples])

#climate data
tsa <- ncvar_get(c.hist, "TSA", start=c(1,1,(syr-1975)*12+1), count=c(nlon,nlat,nyr*12))-273.15 #to celsius
rh  <- ncvar_get(c.hist, "RH2M", start=c(1,1,(syr-1975)*12+1), count=c(nlon,nlat,nyr*12)) #2m relative humidity
ppt <- ncvar_get(c.hist, "RAIN", start=c(1,1,(syr-1975)*12+1), count=c(nlon,nlat,nyr*12))+
       ncvar_get(c.hist, "SNOW", start=c(1,1,(syr-1975)*12+1), count=c(nlon,nlat,nyr*12))
ppt <- ppt*30*3600*24 #mm/s to mm/month
radd <- ncvar_get(c.hist, "FSDSVD", start=c(1,1,(syr-1975)*12+1), count=c(nlon,nlat,nyr*12)) #only vis band is absorbed as PAR
radi <- ncvar_get(c.hist, "FSDSVI", start=c(1,1,(syr-1975)*12+1), count=c(nlon,nlat,nyr*12)) #do not use FSDSND and FSDSNI


samples.m <- samples[,,rep(1:nyr, each=12)] #annual data to monthly
#annual average climate data masked by all crop area
hist.c <- data.table(aod.vis=c((aod.vis2/(dayfoc+1e-10))[samples.m]), aod.volc=c((aod.volc/(dayfoc+1e-10))[samples.m]),
                     ppt=c(ppt[samples.m]), tsa=c(tsa[samples.m]), #2m air temperature
                     rh=c(rh[samples.m]), radd=c(radd[samples.m]), radi=c(radi[samples.m]),
                     Year=rep(rep(syr:fyr,each=12),each=288*192)[samples.m], Month=rep(rep(1:12,time=nyr),each=288*192)[samples.m],
                     Sample=rep(gridmx, time=nyr*12)[samples.m])

clim.yr <- hist.c[,.(aod.vis=mean(aod.vis),aod.volc=mean(aod.volc),tsa=mean(tsa),ppt=mean(ppt),rh=mean(rh),radd=mean(radd),radi=mean(radi)),
                      #tsa.dm=mean(tsa.dm),ppt.dm=mean(ppt.dm),rh.dm=mean(rh.dm),radd.dm=mean(radd.dm),radi.dm=mean(radi.dm)), 
                  by=.(Year,Sample)] #annual average of climate
clim.m <- clim.yr[,.(tsa.m=mean(tsa),ppt.m=mean(ppt),rh.m=mean(rh),radd.m=mean(radd),radi.m=mean(radi)), 
                 by=.(Sample)] #mean per grid sample
clim.yr <- clim.yr[clim.m, .(Year, tsa, ppt, rh, radd, radi, aod.vis, aod.volc,
                   tsa.dm=tsa-i.tsa.m,ppt.dm=ppt-i.ppt.m,rh.dm=rh-i.rh.m,radd.dm=radd-i.radd.m,radi.dm=radi-i.radi.m,
                   radd.pct=(radd-i.radd.m)/i.radd.m*100,radi.pct=(radi-i.radi.m)/i.radi.m*100, 
                   radd.log=log(radd)-log(i.radd.m), radi.log=log(radi)-log(i.radi.m)),
                   by=.EACHI, on=.(Sample)] #de-mean to get the difference climate

library(pracma)
#detrend temperature to remove GHG warming effect
clim.yr <- clim.yr[, tsa.dtr:=detrend(tsa, 'linear'), by=.(Sample)]

df.all <- hist.y[clim.yr, .(area, aod.vis=i.aod.vis, aod.volc=i.aod.volc, 
                            tsa.dtr=i.tsa.dtr, tsa=i.tsa, ppt=i.ppt, rh=i.rh, rg=i.radd+i.radi, radd=i.radd, radi=i.radi,
                            tsa.dif=i.tsa.dm, ppt.dif=i.ppt.dm, rh.dif=i.rh.dm, rg.dif=i.radd.dm+i.radi.dm, 
                            radd.dif=i.radd.dm, radi.dif=i.radi.dm, radd.pct=i.radd.pct, radi.pct=i.radi.pct,
                            radd.log=i.radd.log, radi.log=i.radi.log),
                by=.EACHI, on=.(Year,Sample)]  

df.gl  <- df.all[,.(aod.vis.gl=weighted.mean(aod.vis, area, na.rm = T), aod.volc.gl=weighted.mean(aod.volc, area, na.rm = T),
                    tsa.gl=weighted.mean(tsa, area, na.rm = T), tsa.dtr=weighted.mean(tsa.dtr, area, na.rm = T), 
                    ppt.gl=weighted.mean(ppt, area, na.rm = T),rh.gl=weighted.mean(rh, area, na.rm = T), 
                    rg.gl=weighted.mean(rg, area, na.rm = T),radd.gl=weighted.mean(radd, area, na.rm = T),radi.gl=weighted.mean(radi, area, na.rm = T),
                    tsa.dif=weighted.mean(tsa.dif, area, na.rm = T),ppt.dif=weighted.mean(ppt.dif, area, na.rm = T),rh.dif=weighted.mean(rh.dif, area, na.rm = T),
                    radd.dif=weighted.mean(radd.dif, area, na.rm = T),radi.dif=weighted.mean(radi.dif, area, na.rm = T), #abosulte change in W/m2 
                    radd.pct=weighted.mean(radd.pct, area, na.rm = T),radi.pct=weighted.mean(radi.pct, area, na.rm = T), #% change
                    radd.log=weighted.mean(radd.log, area, na.rm = T),radi.log=weighted.mean(radi.log, area, na.rm = T), #log change
                    area.gl=sum(area, na.rm = T)), by=.(Year)]

#note: need to detrend temperature to remove GHG warming effect!
plot(df.gl$tsa.dif~df.gl$Year,type="l") #Pinatubo eruption, the reduction in insolation led to cooling of about 0.5 °C
plot(df.gl$tsa.dtr~df.gl$Year,type="l")
df.gl[,.(Year, tsa.dif=tsa.dtr-mean(tsa.dtr))] #detrended T change (-0.28°C in 1983, -0.39°C in 1992)
df.gl[Year==1992, tsa.dtr] #-0.39 (<-0.5°C in Proctor et al)
df.gl[Year==1983, tsa.dtr] #-0.28 (crop weighted are larger than global average)
#===total aod is similar to Proctor et al, but volcanic only aod is smaller!
df.gl[Year==1992, aod.vis.gl] #0.295 (total AOD)
df.gl[Year==1983, aod.vis.gl] #0.255 (crop weighted are larger than global average)
df.gl[Year==1992, aod.volc.gl] #0.096 (<0.15 in Proctor et al.)
df.gl[Year==1983, aod.volc.gl] #0.06 (Chichón < Pinatubo)
plot(df.gl$aod.vis.gl~df.gl$Year, type="l") #AOD peaks after 1982 and 1991 very similar to Protocor, but values in other years are non-zero
plot(df.gl$aod.volc.gl~df.gl$Year, type="l") #volcanic only AOD is much smaller than 0.15!!! (Proctor et al exaggerated volcanic impact on crops because they did not remove AOD caused by other aerosols)
#make nicer plots 

#radiation change
plot(df.gl$radd.gl~df.gl$Year,type="l")   #original value global average
plot(df.gl$radi.gl~df.gl$Year,type="l") 
plot(df.gl$radd.dif~df.gl$Year,type="l")  #abosulte decrease -3.5 W/m2 for direct light after Pinatubo in 1992
plot(df.gl$radi.dif~df.gl$Year,type="l")  #diffuse light increase 2.5 W/m2 after Pinatubo in 1992
plot(df.gl$radd.pct~df.gl$Year,type="l")  #direct light only decrease -7.2% after Pinatubo in 1992
plot(df.gl$radi.pct~df.gl$Year,type="l")  #diffuse light only increase +6% after Pinatubo in 1992
df.gl[Year==1992, radd.log]   #direct light decrease -0.063 (log) after Pinatubo in 1992
df.gl[Year==1992, radi.log]   #diffuse light increase +0.058 (log) after Pinatubo in 1992


rm(fsds,solin)
#rm(aod.vis2,aod.volc,tsa,ppt,radd,radi,rh)
nc_close(c.hist)
nc_close(y.hist)




#============Figures for response to reviewers (R1)

#==========draw all climate data in one panel
#svg(paste0(wd1,"/crop_figures/revision/Fig.R1.Historical_AOD_climate.svg"), width =6, height = 8) 
jpeg(paste0(wd1,"/figures/supplementary/Fig.R1.Historical_AOD_climate_global.jpeg"), units = "in", width =6, height=10, res = 600)

par(mar=c(1, 5, 1, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
split.screen(c(6,1))
screen(1) # prepare screen 1 for output
plot(aod.vis~Year, data=solar.gl, xlim=c(1979,2009), type="l", bty="n", xaxt="n", xlab="", ylab="NorESM1-ME TAUE550 \n(550nm)") #AOD peaks after 1982 and 1991 very similar to Protocor, but values in other years are non-zero
abline(h=mean(solar.gl[(Year!=1992 & Year!=1983)]$aod.vis), lty="longdash") #reference: mean excluding 1983 and 1992

screen(2) # prepare screen 1 for output
par(mar=c(1, 5, 0, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
plot(aod.vis2~Year, data=solar.gl2, type="l", bty="n", xaxt="n", xlab="", ylab="NorESM2 AOD_VIS \n(0.442-0.625um)") #AOD peaks after 1982 and 1991 very similar to Protocor, but values in other years are non-zero
abline(h=mean(solar.gl2[(Year!=1992 & Year!=1983)]$aod.vis2), lty="longdash") #reference: mean excluding 1983 and 1992

screen(n =3, new = TRUE)
par(mar=c(0, 5, 0, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
plot(aod.volc~Year, data=solar.gl2, type="l", bty="n", xaxt="n", yaxt="n", xlab="", ylab="CMIP6 volcanic-only \nAOD (0.442-0.625um)",
     ylim=c(0,0.1)) 
abline(h=mean(solar.gl2[(Year!=1992 & Year!=1983)]$aod.volc), lty="longdash") #reference: mean excluding 1983 and 1992
myTicks = axTicks(2)
axis(2, at = myTicks, labels = formatC(myTicks, format = 'g', drop0trailing = TRUE))

screen(n =4)
par(mar=c(0, 5, 0, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
plot(radd.pct~Year, data=df.gl, type="l", bty="n", xaxt="n", xlab="", ylab="Direct radiation \ndeviation (%)",
     ylim=c(-7,5)) 
abline(h=mean(df.gl$radd.pct), lty="longdash")

screen(n =5)
par(mar=c(0, 5, 0, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
plot(radi.pct~Year, data=df.gl, type="l", bty="n", xaxt="n", xlab="", ylab="Diffuse radiation \ndeviation (%)",
     ylim=c(-5,6)) 
abline(h=mean(df.gl$radi.pct), lty="longdash")

screen(n =6)
par(mar=c(2.5, 5, 0, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9) #mpg: margins for axis title, label and axis line
plot(tsa.dtr~Year, data=df.gl, type="l", bty="n", xaxt="n", yaxt="n", xlab="Year", ylab="Global average \ntemperature (°C)",
     ylim=c(-0.6,0.6)) 
abline(h=mean(df.gl$tsa.dtr), lty="longdash") #reference: mean excluding 1983 and 1992
myTicks = axTicks(2)
axis(2, at = myTicks, labels = formatC(myTicks, format = 'g', drop0trailing = TRUE))

close.screen(all = TRUE)

par(mar=c(2.3, 4.6, 0.5, 1), mgp=c(1.4,0.5,0), las=1, new=TRUE) #overlay plot to the above split plot
plot(tsa.dtr~Year, data=df.gl, type="n", bty="n", xaxt="l", yaxt="n", xlab="Year", ylab="", cex=0.9)
abline(v=1982, lty="dotted") 
#text(x=1982, y=-0.1, "El Chichón", pos=2, srt=90) #mark the real year when volcanos happened
text(x=1982, y=-0.2, "El Chichón", pos=2, srt=90) #mark the real year when volcanos happened
abline(v=1991, lty="dotted") 
#text(x=1991, y=-0.1, "Pinatubo", pos=2, srt=90)
text(x=1991, y=-0.2, "Pinatubo", pos=2, srt=90)

dev.off()



#========global average (crop area weighted)
#svg(paste0(wd1,"/crop_figures/revision/Fig.R1.Historical_AOD_climate.svg"), width =6, height = 8) 
jpeg(paste0(wd1,"/figures/supplementary/Fig.R1.Historical_AOD_climate_croparea.jpeg"), units = "in", width =6, height=10, res = 600)

par(mar=c(1, 5, 1, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
split.screen(c(5,1))
screen(1) # prepare screen 1 for output
par(mar=c(0, 5, 0, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
plot(aod.vis.gl~Year, data=df.gl, type="l", bty="n", xaxt="n", xlab="", ylab="Atmospheric AOD \n(0.442-0.625um)") #AOD peaks after 1982 and 1991 very similar to Protocor, but values in other years are non-zero
abline(h=mean(df.gl[(Year!=1992 & Year!=1983)]$aod.vis.gl), lty="longdash") #reference: mean excluding 1983 and 1992

screen(n =2, new = TRUE)
par(mar=c(0, 5, 0, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
plot(aod.volc.gl~Year, data=df.gl, type="l", bty="n", xaxt="n", yaxt="n", xlab="", ylab="CMIP6 volcanic-only \nAOD (0.442-0.625um)",
     ylim=c(0,0.1)) 
abline(h=mean(df.gl[(Year!=1992 & Year!=1983)]$aod.volc.gl), lty="longdash") #reference: mean excluding 1983 and 1992
myTicks = axTicks(2)
axis(2, at = myTicks, labels = formatC(myTicks, format = 'g', drop0trailing = TRUE))

screen(n =3)
par(mar=c(0, 5, 0, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
plot(radd.pct~Year, data=df.gl, type="l", bty="n", xaxt="n", xlab="", ylab="Direct radiation \ndeviation (%)",
     ylim=c(-7,5)) 
abline(h=mean(df.gl$radd.pct), lty="longdash")

screen(n =4)
par(mar=c(0, 5, 0, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9)
plot(radi.pct~Year, data=df.gl, type="l", bty="n", xaxt="n", xlab="", ylab="Diffuse radiation \ndeviation (%)",
     ylim=c(-5,6)) 
abline(h=mean(df.gl$radi.pct), lty="longdash")

screen(n =5)
par(mar=c(2.5, 5, 0, 1) + 0.1, mgp=c(2.5,0.5,0), las=1, cex =0.9) #mpg: margins for axis title, label and axis line
plot(tsa.dtr~Year, data=df.gl, type="l", bty="n", xaxt="n", yaxt="n", xlab="Year", ylab="Global average \ntemperature (°C)",
     ylim=c(-0.6,0.6)) 
abline(h=mean(df.gl$tsa.dtr), lty="longdash") #reference: mean excluding 1983 and 1992
myTicks = axTicks(2)
axis(2, at = myTicks, labels = formatC(myTicks, format = 'g', drop0trailing = TRUE))

close.screen(all = TRUE)

par(mar=c(2.3, 4.6, 0.5, 1), mgp=c(1.4,0.5,0), las=1, new=TRUE) #overlay plot to the above split plot
plot(tsa.dtr~Year, data=df.gl, type="n", bty="n", xaxt="l", yaxt="n", xlab="Year", ylab="", cex=0.9)
abline(v=1982, lty="dotted") 
text(x=1982, y=-0.1, "El Chichón", pos=2, srt=90) #mark the real year when volcanos happened
abline(v=1991, lty="dotted") 
text(x=1991, y=-0.1, "Pinatubo", pos=2, srt=90)

dev.off()





#=======growing season climate data masked by crop-specific cultivation area
crops = c("Corn", "Wheat", "Soybean", "Cotton", "Rice", "Sugarcane", "Tropical Corn", "Tropical Soybean")
lon <- ncvar_get(c.hist,"lon")
lat <- ncvar_get(c.hist,"lat")
nlon <- length(lon); nlat<-length(lat)
latmx <- matrix(rep(lat,each=nlon),nrow=nlon,ncol=nlat) #matrix of latitude
gridmx <- array(1:(nlat*nlon),list(nlon,nlat)) #matrix of grid index

df.hist.Gr<-data.table()
for (n in 1:8) {
  for (i in 0:1){
    if (i==0){
      ncft = seq(1,15,2)[n]  #rainfed
      irrigated=FALSE
    }else{
      ncft = seq(2,16,2)[n]  #irrigated
      irrigated=TRUE
    }  
    if (crops[n]=="Rice") {
      GrowNH=c(0,1,1,1,1,1,1,0,0,0,0,0) #Rice plants earlier, growing from Feb-Jun
      GrowSH=c(1,0,0,0,0,0,0,1,1,1,1,1)
    } else if (crops[n]=="Sugarcane"){ #sugarcane grows 10 months, planting and harvest dates similar, so use annual climate
      GrowNH=c(1,1,1,1,1,1,1,1,1,1,1,1)
      GrowSH=c(1,1,1,1,1,1,1,1,1,1,1,1)
    } else {
      GrowNH=c(0,0,0,0,1,1,1,1,1,1,0,0) #most crops grow 150 days from April/June to Aug/Oct
      GrowSH=c(1,1,1,1,0,0,0,0,0,0,1,1) #use average May-Oct for NH, Nov-Apr for SH, as in Levis et al. 2016
    }
    
    area <- cft.area[,,ncft,]
    #all grid-year samples with >1000ha of a cft or 0.1% of grid cell in RCP85
    samples <- area>=0.001 #Mha -> 1000ha
    samples[is.na(samples)] <- FALSE
    #samplesize <- apply(samples, 3, function(x)sum(array(1, list(nlon,nlat))[x], na.rm = T)) #different number of samples per year
    #samplesize #around 8100~8200 grid cells with crops
    
    samples.m <- samples[,,rep(1:nyr, each=12)] #annual data to monthly
    area.m <- area[,,rep(1:nyr, each=12)] #annual data to monthly
    hist.c <- data.table(aod.vis=c(aod.vis2[samples.m]), aod.volc=c(aod.volc[samples.m]),
                         ppt=c(ppt[samples.m]), tsa=c(tsa[samples.m]), #2m air temperature
                         rh=c(rh[samples.m]), radd=c(radd[samples.m]), radi=c(radi[samples.m]),
                         Year=rep(rep(syr:fyr,each=12),each=288*192)[samples.m], Month=rep(rep(1:12,time=nyr),each=288*192)[samples.m],
                         GrowingN=rep(rep(GrowNH,time=nyr),each=288*192)[samples.m], #growing season filter: April-Oct (use LAI>0.5)
                         GrowingS=rep(rep(GrowSH,time=nyr),each=288*192)[samples.m], #growing season filter: (most crops have only 5 months growing season)
                         area = c(area.m[samples.m]), #crop specific area
                         Sample=rep(gridmx, time=nyr*12)[samples.m], Lat=rep(latmx, time=nyr*12)[samples.m])
    
    #===Option 1: calculate growing season climate for ech crop according to their planting date and max harvest date
    hist.NH <- hist.c[Lat>=0,.(tsa=mean(tsa), aod.vis=mean(aod.vis), aod.volc=mean(aod.volc), 
                               rh=mean(rh), ppt=mean(ppt),radd=mean(radd),radi=mean(radi),
                               area=mean(area)), by=.(Sample,Year,GrowingN)] #growing season mean using NH filter
    hist.SH <- hist.c[Lat<0,.(tsa=mean(tsa), aod.vis=mean(aod.vis), aod.volc=mean(aod.volc),
                              rh=mean(rh), ppt=mean(ppt),radd=mean(radd),radi=mean(radi), 
                              area=mean(area)), by=.(Sample,Year,GrowingS)] #SH filter
    hist.G1 <- rbind(hist.NH[GrowingN==1], hist.SH[GrowingS==1], fill=T)
    hist.G1 <- hist.G1[order(Sample,Year),]
    #hist.G1 <- hist.G1[hist.y, on=.(Sample,Year)] #merge yield and climate data
    hist.G1$Crop=crops[n]
    hist.G1$Irrig=irrigated
    
    ##=====merge output data of all crop types (8x2) 
    df.hist.Gr <- rbind(df.hist.Gr, hist.G1, fill=TRUE)
    rm(hist.G1, hist.NH, hist.SH)
  }}

#detrend temperature to remove GHG warming effect
df.hist.Gr <- df.hist.Gr[, tsa.dtr:=detrend(tsa, 'linear'), by=.(Sample)] #detrended T per crop will remove too many grid cells
df.hist.m <- df.hist.Gr[,.(tsa.m=mean(tsa), tsa.dtr.m=mean(tsa.dtr,na.rm=TRUE), 
                           ppt.m=mean(ppt),rh.m=mean(rh),radd.m=mean(radd),radi.m=mean(radi)), 
                  by=.(Sample,Crop,Irrig)] #mean per grid sample per crop

df.hist.Gr<- df.hist.Gr[df.hist.m, .(Year, area, tsa, tsa.dtr, ppt, rh, radd, radi, aod.vis, aod.volc,
                             tsa.dif=tsa-i.tsa.m, tsa.dif1=tsa.dtr-i.tsa.dtr.m,
                             ppt.dif=ppt-i.ppt.m,rh.dif=rh-i.rh.m,radd.dif=radd-i.radd.m,radi.dif=radi-i.radi.m,
                             radd.pct=(radd-i.radd.m)/i.radd.m*100,radi.pct=(radi-i.radi.m)/i.radi.m*100, 
                             radd.log=log(radd)-log(i.radd.m), radi.log=log(radi)-log(i.radi.m)),
                   by=.EACHI, on=.(Sample,Crop,Irrig)] #de-mean to get the difference climate

rm(cft.area,area.m,area)
rm(aod.vis2,aod.volc,tsa,ppt,radd,radi,rh)
rm(samples, samples.m)


#================answers=================
#Volcano-only AOD is also much lower in the model (+0.048 for 1992, +0.03 for 1983) than in Proctor et al. (+0.15)
#and volcanic impact on direct and diffuse radiation is much smaller in NorESM2 than in Proctor et al. (direct -21%, diffuse +20% for Pinatubo)
#By the end of 21st century under SAI, changes in direct (-13%) and diffuse (+23%) radiation are non-symmetrical, thus resulting in little net insolation effect 


#Reviewer#3: could an estimate of the yield response in CLM attributable to such radiative changes be made using the results of the sensitivity analysis?
#option 1: apply the coefficients of our sensitivity analysis to the radiative differences between modeled SAI in our study and volcanic data in Proctor et al. to estimate the difference in yield responses between the two studies
#option 2: apply the coefficients of our sensitivity analysis to the above radiative changes from Proctor et al. (assumed uniformly across grid cells) to estimate yield response 


#1. global estimate with mean coefficients
load(paste0(wd1,"/data/RData/MLR-pergridregression.Rdata"))
coeff.i.sai.m <- coeff.i.samples[Scenario=="SAI - RCP85",.(Intercept=mean(Intercept), T=mean(T), RD=mean(RD), RI=mean(RI)), by=.(Scenario,Crop)]
coeff.r.sai.m <- coeff.r.samples[Scenario=="SAI - RCP85",.(Intercept=mean(Intercept), T=mean(T), RD=mean(RD), RI=mean(RI), P=mean(P), RH=mean(RH) ), by=.(Scenario,Crop)]

Pinatubo.clim <- df.gl[Year==1992,]
Chichón.clim <- df.gl[Year==1983,]
corn.r <- coeff.r.sai.m[Crop=="Corn"]
Pinatubo.log <- corn.r$Intercept+corn.r$T*Pinatubo.clim$tsa.dif+corn.r$RD*Pinatubo.clim$radd.dif+corn.r$RI*Pinatubo.clim$radi.dif + 
                corn.r$P*Pinatubo.clim$ppt.dif + corn.r$RH*Pinatubo.clim$rh.dif
Pinatubo.eff <- (exp(Pinatubo.log)-1)*100
(T.eff <- (exp(corn.r$T*Pinatubo.clim$tsa.dif)-1)*100) #+0.15%, T only decreased -0.236°C in 1992 (Proctor says global cooling of about 0.5 °C)
(T.eff <- (exp(corn.r$T*(-0.5))-1)*100) #+0.45%, very small temperature effect by Pinatubo even if temperature decrease by -0.5 °C
(R.eff <- (exp(corn.r$RD*Pinatubo.clim$radd.dif+corn.r$RI*Pinatubo.clim$radi.dif)-1)*100) #+1%, positive insolation effect
(R.eff <- (exp(corn.r$RD*Pinatubo.clim$radd.dif*3+corn.r$RI*Pinatubo.clim$radi.dif*3)-1)*100) #scaled up by 3 times to match radiative change of Pinatubo

Chichón.log <- corn.r$Intercept+corn.r$T*Chichón.clim$tsa.dif+corn.r$RD*Chichón.clim$radd.dif+corn.r$RI*Chichón.clim$radi.dif + 
               corn.r$P*Chichón.clim$ppt.dif + corn.r$RH*Chichón.clim$rh.dif
Chichón.eff <- (exp(Chichón.log)-1)*100
(T.eff <- (exp(corn.r$T*Chichón.clim$tsa.dif)-1)*100) #0.84%; 
#stronger cooling (-0.91°C) after Chichón with reference to the mean of 1980-2014 (need to detrend T to remove global warming impact)
#detrended T change (-0.188°C in 1983, -0.113°C in 1992)
(R.eff <- (exp(corn.r$RD*Chichón.clim$radd.dif+corn.r$RI*Chichón.clim$radi.dif)-1)*100) #1.6%


#2. estimate with per-grid-cell coefficients for each crop type (here take corn for example, because df.all does not have crop-specific climate data)
coeff.i.sai <- coeff.i.samples[Scenario=="SAI - RCP85" & Crop=="Corn",]
coeff.r.sai <- coeff.r.samples[Scenario=="SAI - RCP85" & Crop=="Corn",]
unique(coeff.i.sai$Crop)

effect.i.sai <- coeff.i.sai[df.all, .(Intercept, T.log=T*i.tsa.dif, RD.log=RD*i.radd.dif, RI.log=RI*i.radi.dif, 
                                      R.log=RD*i.radd.dif+RI*i.radi.dif, P.log=0, RH.log=0, 
                                      area=i.area, Year=i.Year), by=.EACHI, on=.(Crop,Sample)]
effect.r.sai <- coeff.r.sai[df.all, .(Intercept, T.log=T*i.tsa.dif, RD.log=RD*i.radd.dif, RI.log=RI*i.radi.dif, 
                                      R.log=RD*i.radd.dif+RI*i.radi.dif, RH.log=RH*i.rh.dif,P.log=P*i.ppt.dif,
                                      area=i.area, Year=i.Year), by=.EACHI, on=.(Crop,Sample)]
effect.i.sai <- effect.i.sai[,Tot.log:=T.log+RD.log+RI.log+Intercept] 
effect.r.sai <- effect.r.sai[,Tot.log:=T.log+RD.log+RI.log+P.log+RH.log+Intercept]
effect.i.sai <- effect.i.sai[!is.na(Tot.log)]
effect.r.sai <- effect.r.sai[!is.na(Tot.log)]
effect.i.sai$Irrig <- TRUE
effect.r.sai$Irrig <- FALSE
effect.sai <- rbind(effect.i.sai, effect.r.sai)


eff.r.wm <- effect.r.sai[, .(T.log=weighted.mean(T.log, area, na.rm = T), RD.log=weighted.mean(RD.log, area, na.rm = T),RI.log=weighted.mean(RI.log, area, na.rm = T),
                             R.log=weighted.mean(R.log, area, na.rm = T), P.log=weighted.mean(P.log, area, na.rm = T),RH.log=weighted.mean(RH.log, area, na.rm = T),
                           Tot.log=weighted.mean(Tot.log, area, na.rm = T), #Y.mod.dif=weighted.mean(Y.mod.dif, area, na.rm = T), 
                           #Res.log=weighted.mean(Res.log, area, na.rm = T),
                           area=sum(area, na.rm=T)),
                       by =.(Year)]

eff.r.pct <- eff.r.wm[,.(T.eff=(exp(T.log)-1)*100, R.eff=(exp(R.log)-1)*100, RD.eff=(exp(RD.log)-1)*100, RI.eff=(exp(RI.log)-1)*100, 
                         P.eff=(exp(P.log)-1)*100, RH.eff=(exp(RH.log)-1)*100,
                         Tot.eff=(exp(Tot.log)-1)*100), by =.(Year)]
eff.r.pct[Year==1992] #+0.49%
eff.r.pct[Year==1983] #+1.14% total radiation effect of Chichón indeed higher than Pinatubo
#total climate effect on maize: 2.5% for Pinatubo, 2.06% for Chichón (similar to the mean value in Fig.4e in Proctor et al.)


#=============next step=============== 
#R effect slight positive because radiative changes in NorESM (-7.2% direct, +6% diffuse in 1992) are smaller than Pinatubo eruption (-21% direct, +20% diffuse) but both are more symmetrical than future SAI
#T effect too small for 1992 and 1983 because cooling by these historical volcanos were much smaller than future SAI experiments
#detrended global T change (-0.188°C in 1983, -0.113°C in 1992)
#==more realistic way to estimate volcanic effect: scale up radiative changes of SAI experiment to match % change as of Pinatubo
#==and then apply SAI coefficients to the new radiative changes for the past or future

coeff.r.samples.m[Scenario=="SAI - RCP85", RI/RD, by=.(Scenario,Crop)] 
coeff.r.samples.m[Scenario=="SAI - RCP85", mean(RI/RD), by=.(Scenario)] #average 1.19 across crop types
#In our SAI crops diffuse:direct effect ratio, range from 0.54 of wheat to 1.76 of corn, average 1.19, very similar to RI < 1.2RD (g2 < 1.2g1)
#Section III.5 Proctor et al:
#"Our finding that the insolation effect is negative suggests that g1 > g2w. Calculation of wss using our insolation data, suggests that w ~ 0.85. 
#"This means that g2 < 1.2g1, meaning that diffuse light is less than 1.2 times more efficiently used for the production of yield than direct light."


#==example 1 (assume the % changes in direct (-21%) and diffuse (+20%) light as in Proctor et al.)
#==apply our SAI coefficients to rescaled radiative changes during 1979–2009, while other climate variables use HIST
load(paste0(wd1,"/data/RData/Fig.S3.MLR-pergrid-rad-ppt.log.Rdata"))
coeff.i.sai <- coeff.i.samples3[Scenario=="SAI - RCP85"] #use log change as predictor for RADD, RADI
coeff.r.sai <- coeff.r.samples3[Scenario=="SAI - RCP85"]

str(df.hist.Gr) #rescale df.hist.Gr so that the radd.log and radi.log become similar to -21% direct, +20% diffuse
df.hist.Gr[Year==1992, weighted.mean(tsa.dif,area,na.rm=T)] #orginal T change, -0.48 degree (close to -0.5 in Proctor)
df.hist.Gr[Year==1992, weighted.mean(tsa.dif1,area,na.rm=T)]  #detrended T change, -0.45 (too many grid cells removed due to small sample size for detrending)
(radd.log.m <- df.hist.Gr[Year==1992, weighted.mean(radd.log,area,na.rm=T)]) #direct light decrease -0.058 (-5.7%) after Pinatubo in 1992
(radi.log.m <- df.hist.Gr[Year==1992, weighted.mean(radi.log,area,na.rm=T)]) #diffuse light increase +0.06 (+6.3%) after Pinatubo in 1992

#1) rescale radd/radi change by log(1-0.21)/-0.058 =4., log(1+0.2)/0.06=3. to match historical radiative change (-21% direct, +20% diffuse after Pinatubo) in Proctor et al.
(scaler1 <- log(1-21/100)/radd.log.m ); 
(scaler2 <- log(1+20/100)/radi.log.m )
effect.i.hist <- coeff.i.sai[df.hist.Gr[Irrig==TRUE], 
                             .(Intercept, T.log=T*i.tsa.dif, P.log=0, RH.log=0, 
                               #RD.log=RD*i.radd.dif, RI.log=RI*i.radi.dif, R.log=RD*i.radd.dif+RI*i.radi.dif, 
                               #RD.log=RD*i.radd.log, RI.log=RI*i.radi.log, R.log=RD*i.radd.log+RI*i.radi.log, #when coefficients from log RADD/RADI as predictors
                               RD.log=RD*i.radd.log*scaler1, RI.log=RI*i.radi.log*scaler2, R.log=RD*i.radd.log*scaler1+RI*i.radi.log*scaler2, #rescaled
                               area=i.area, Year=i.Year), by=.EACHI, on=.(Crop,Sample)]
effect.r.hist <- coeff.r.sai[df.hist.Gr[Irrig==FALSE], 
                             .(Intercept, T.log=T*i.tsa.dif, RH.log=RH*i.rh.dif,P.log=P*i.ppt.dif,
                               #RD.log=RD*i.radd.dif, RI.log=RI*i.radi.dif, R.log=RD*i.radd.dif+RI*i.radi.dif,
                               #RD.log=RD*i.radd.log, RI.log=RI*i.radi.log, R.log=RD*i.radd.log+RI*i.radi.log,
                               RD.log=RD*i.radd.log*scaler1, RI.log=RI*i.radi.log*scaler2, R.log=RD*i.radd.log*scaler1+RI*i.radi.log*scaler2, #rescaled
                               area=i.area, Year=i.Year), by=.EACHI, on=.(Crop,Sample)]
effect.i.hist <- effect.i.hist[,Tot.log:=T.log+RD.log+RI.log+Intercept] 
effect.r.hist <- effect.r.hist[,Tot.log:=T.log+RD.log+RI.log+P.log+RH.log+Intercept]
effect.i.hist <- effect.i.hist[!is.na(Tot.log)]
effect.r.hist <- effect.r.hist[!is.na(Tot.log)]
effect.i.hist$Irrig <- TRUE
effect.r.hist$Irrig <- FALSE
effect.hist <- rbind(effect.i.hist, effect.r.hist)

eff.hist.wm <- effect.hist[, .(T.log=weighted.mean(T.log, area, na.rm = T), R.log=weighted.mean(R.log, area, na.rm = T),
                                 RD.log=weighted.mean(RD.log, area, na.rm = T),RI.log=weighted.mean(RI.log, area, na.rm = T),
                                 RH.log=weighted.mean(RH.log, area, na.rm = T), P.log=weighted.mean(P.log, area, na.rm = T), 
                                 Tot.log=weighted.mean(Tot.log, area, na.rm = T), #Y.mod.log=weighted.mean(Y.mod.log, area, na.rm = T), 
                                 area=sum(area, na.rm=T)),
                             by =.(Crop,Irrig,Year)]
unique(eff.hist.wm$Year);unique(eff.hist.wm$Crop)

#merge effects of 8 crops to 6 crops and merge irrig/rainfed
eff.hist.wm[Crop=="Tropical Corn", Crop:="Corn"]
eff.hist.wm[Crop=="Tropical Soybean", Crop:="Soybean"]

#per crop effect: merge irrig/rainfed, weight by area only (no data for base yield of each crop type)
eff.hist <- eff.hist.wm[, .(T.log=weighted.mean(T.log, area, na.rm = T), R.log=weighted.mean(R.log, area, na.rm = T),
                             RD.log=weighted.mean(RD.log, area, na.rm = T),RI.log=weighted.mean(RI.log, area, na.rm = T),
                             P.log=weighted.mean(P.log, area, na.rm = T), RH.log=weighted.mean(RH.log, area, na.rm = T),
                             Tot.log=weighted.mean(Tot.log, area, na.rm = T), 
                             area=sum(area, na.rm = T)), by =.(Crop,Year)]
#translate log to % effect
eff.hist.pct <- eff.hist[, .( T.eff=(exp(T.log)-1)*100, R.eff=(exp(R.log)-1)*100, 
                              RD.eff=(exp(RD.log)-1)*100, RI.eff=(exp(RI.log)-1)*100, P.eff=(exp(P.log)-1)*100, RH.eff=(exp(RH.log)-1)*100,
                              Tot.eff=(exp(Tot.log)-1)*100, area), by=.(Crop,Year)]
eff.hist.pct[Year==1983]
eff.hist.pct[Year==1992] #negative insolation effect for most crops when log changes in RADD, RADI are used as predictors
#net insolation effect (Corn -1.0%, Rice -0.53%, Soy 0.36%, Wheat -0.84%), direct and diffuse effects largely cancel out each other
#T.eff (Corn +3.9%, Rice 1.5%, Soy 9.9%, Wheat -1.4%) with T changed -0.48°C in 1992 from the mean of 1979-2009
#Tot.eff (Corn +3.7%, Rice -1.3%, Soy 5.5%, Wheat 7.5%) in 1992, Rice (-5%) and Wheat (-6.2%) negative in 1983

#===after rescaling radd.log and radi.log to match radiative change in Proctor (-21% direct, +20% diffuse)
#net insolation effect become more negative for Pinatubo (Corn -8.9%, Rice -6.7%, Soy -3.6%, Wheat -2.7%), very similar to Fig.3b in Proctor et al.
#net insolation effect for El Chichón more positive (Corn -1.4%, Rice 5.7%, Soy 15.2%, Wheat 0.9%)
#Tot.eff (Corn -4.6%, Rice -7.5%, Soy 1.4%, Wheat 5.5%) in 1992, (Corn 1.3%, Rice -1.4%, Soy 15.3%, Wheat -5.2%) in 1983






#2) rescale our SAI - RCP85 radd/radi changes to match future radiative changes (-11% direct, +10% diffuse) of 2050-2069 with SAI treatment in Proctor et al.
dAOD = 0.084 # for the SAI - RCP45 stablizing experiment in Proctor et al.
(exp(-1.406*dAOD)-1)*100; (exp(1.171*dAOD)-1)*100 #-11% () direct / +10.3% diffuse if using AOD coefficient from Pinatubo 
#(exp(-1.039*dAOD)-1)*100; (exp(2.063*dAOD)-1)*100 #-8% direct / +19% diffuse if using AOD coefficient from El Chichón 
#these radiative changes (-11% direct +10.3 diffuse) are very similar to unit changes in the below figure
#so our SAI coefficients could reproduce the negative insolation effect seen in Proctor et al. if direct/diffuse changes are symmetrical
dif.unit <- dif.8crop.s[Scenario=="SAI - RCP85",
                        .(tsa.unit=-0.88, radd.unit=log(1-1.406*dAOD), radi.unit=log(1+1.171*dAOD), 
                          yield0, area.wt,Year,Crop,Irrig,Sample)]

eff.i.sai <- coeff.i.sai[dif.unit[Irrig==TRUE],                                    
                          .(T.log=T*i.tsa.unit, RD.log=RD*i.radd.unit, RI.log=RI*i.radi.unit, 
                            R.log=RD*i.radd.unit+RI*i.radi.unit, 
                            yield0=i.yield0, area=i.area.wt, Year=i.Year), by=.EACHI, on=.(Crop,Sample)]
eff.r.sai <- coeff.r.sai[dif.unit[Irrig==FALSE], 
                          .(T.log=T*i.tsa.unit, R.log=RD*i.radd.unit+RI*i.radi.unit, 
                            RD.log=RD*i.radd.unit, RI.log=RI*i.radi.unit, 
                            yield0=i.yield0,area=i.area.wt, Year=i.Year), by=.EACHI, on=.(Crop,Sample)]
eff.i.sai$Irrig <- TRUE
eff.r.sai$Irrig <- FALSE
eff.sai <- rbind(eff.i.sai, eff.r.sai)

eff.sai.wm <- eff.sai[, .(T.log=weighted.mean(T.log, area, na.rm = T), R.log=weighted.mean(R.log, area, na.rm = T),
                            RD.log=weighted.mean(RD.log, area, na.rm = T),RI.log=weighted.mean(RI.log, area, na.rm = T),
                            yield0=weighted.mean(yield0, area, na.rm = T), area=sum(area, na.rm=T)),
                        by =.(Crop,Irrig,Year)]

eff.sai.pct <- eff.sai.wm[Year>=2050 & Year<=2069, .( T.eff=(exp(T.log)-1)*100, R.eff=(exp(R.log)-1)*100, 
                                   RD.eff=(exp(RD.log)-1)*100, RI.eff=(exp(RI.log)-1)*100, yield0, area), 
                          by=.(Crop,Irrig)]
#merge effects of 8 crops to 6 crops and merge irrig/rainfed
eff.sai.pct[Crop=="Tropical Corn", Crop:="Corn"]
eff.sai.pct[Crop=="Tropical Soybean", Crop:="Soybean"]

#per crop effect: merge irrig/rainfed, weight by area and base yield
eff.sai <- eff.sai.pct[, .(T.eff=weighted.mean(T.eff, yield0*area, na.rm = T), R.eff=weighted.mean(R.eff, yield0*area, na.rm = T),
                             RD.eff=weighted.mean(RD.eff, yield0*area, na.rm = T),RI.eff=weighted.mean(RI.eff, yield0*area, na.rm = T),
                             area=sum(area, na.rm = T)), by =.(Crop)]

eff.sai #net insolation effect negative for all crops (Corn -5.5%, Rice -3.8%, Soy -8.2%, Wheat -3.2%), very similar to Fig.4e in Proctor et al.


















